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Quantum  Simulation  of  Heterogeneous  Electron  Transfer  Free 
Energies  at  the  Water*Metal  Interface 

Jay  B.  Straus*  and  Gregory  A.  Voth 

Department  of  Chemutry,  Univeniiy  of  Penntylvattia,  Philaddphia,  PA  19104~6323,  USA 

Abstract 

The  free  energy  for  electnm  transfer  is  studied  for  the  Fe^'*'/Fe^  charge 
transfer  system  at  the  water>Pt(lll)  interface.  Classical  diabatic  free  en¬ 
ergy  curves  are  calculated  along  with  with  an  adiabatic  curve  baaed  mi  the 
Anderson-Newns  Hamiltonian.  Reactive  flux  calculations  are  then  performed 
on  this  curve  to  determine  the  effect  rtf  recrossings  on  the  classical  rate  con¬ 
stant.  These  effects  are  not  found  to  be  large  (k  ~  0.6).  Uie  solvent  model  is 
then  extended  to  a  quantum  mechanical  path  integral  version  and  the  quan¬ 
tum  aiflabatic  free  energy  curves  are  calculated.  The  resulting  quantum  ef¬ 
fects  are  found  to  bequite  rigniflcant,  illustrating  that  the  same  dectrode 
overpotential  does  not  necessarily  result  in  the  same  free  energy  curves  for 
the  classical  and  quantum  medianical  solvent  models.  These  results  suggest 
that  classical  modds  for  water  mey  not  be  adequate,  or  at  least  need  to  be 
modified,  for  accurate  simulations  of  heterogeneous  dectron  transfer. 
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I.  INTRODUCTION 


EleGtron  tmufer  (ET)  is  a  process  of  fundamental  importance  in  chemistry,  physics,  and 
bi<^ogy  [1,2].  The  field  has  seen  substantial  theoretical  research  activity,  including  computer 
simulations  of  a  variety  of  homogeneous  electron  transfer  systems,  such  as  charge  separation 
in  an  ion  pair  in  a  polar  solvent  [3-7]  as  wdl  as  electron  transfer  in  the  homogeneous  ferrous- 
ferric  qrstem.  [8-10]  Recently,  however,  there  has  been  an  emerging  focus  on  heterogeneous 
ET  sjrstems  [11-13]  firom  the  pmnt  of  view  of  computer  simulation.  It  is  such  interfadal 
charge  transfer  phenomenon  which  are  important,  for  example,  at  a  metal  electrode  in^ 
mersed  in  an  electrolyte  solution.  A  better  understanding  of  heterogeneous  ET  processes 
will  be  relevant  to  many  areas  of  technology,  ranging  from  electrodiemistry  to  ccnrosion  to 
battery  technology.  The  present  authors  have  previously  performed  simulations  examining 
the  firee  energy  curves  for  a  moddi  ion  in  water  near  a  Pt(lll)  surfoce.  [11]  That  work  showed 
an  insensitivity  of  the  solvent  diabatic  free  energy  curves  to  surface-induced  solvent  density 
inhomogeoidties.  It  also  illustrated  a  substantial  adherence  to  Marcus  parabolic  behavior  [1] 
over  the  wide  *^ormal”  free  energy  range,  as  well  as  a  significant  deviation  in  the  inverted 
legone.  The  present  work  significantly  extends  this  area  of  research  with  the  focus  being  on 
a  more  realistic  ferrous-ferric  charge  transfer  system  at  a  platinum  dectrode.  In  additimi  to 
finding  the  diabatic  free  energy  curves,  a  classical  adiabatic  curve  is  computed  using  a  modd 
for  heterogeneous  ET  based  on  the  Anderson-Newns  Hamiltonian.  [14r-17]  Classical  reactive 
flax  calculations  are  then  performed  for  this  system  to  determine  the  effect  of  transition 
state  rectossings  within  the  Marcus  theoretical  framework  for  the  classical  adiabatic  rate 
constant.  This  latter  work  is  related  to  the  work  of  Rose  and  Benjamin  except  that  those 
authors  used  a  sinq>Ier  two-state  description  of  the  heterogeneous  ET  Hamiltonian.  [13] 

In  a  significant  departure  from  previous  simulation  studies  of  heterogeneous  ET,  a  path 
integral  modd  is  developed  in  order  to  probe  the  effects  of  quantization  of  the  water  nudear 
motion.  This  model  is  employed  to  determine  the  quantum  adiabatic  free  energy  curves 
within  the  context  of  path  integral  quantum  tranrition  state  theory.  [18]  The  important 
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effects  of  water  quantisation  on  the  ET  free  energies  is  then  assessed,  with  special  attention 
placed  on  their  magnitude  with  respect  to  the  size  of  the  classical  transition  state  recrossing 
effects.  It  is  found  that  the  effect  of  quantizing  the  water  modes  is  by  far  the  most  important 
one,  leading  to  the  conclusion  that  the  effects  of  water  quantization  must  be  considered  in 
order  for  simulations  of  heterogeneous  ET  to  be  accurate. 

n.  METHODS 

A.  The  Anderson-Newns  Hamiltonian 

Charge  transfer  between  an  ion  in  solution  and  a  metal  electrode  is  represented  in  the 
present  researdi  by  a  version  [19-21]  of  an  Hamiltonian  referred  to  in  the  literature  as  the 
Andersmi'Newns  Hamiltonian  [14-17].  It  can  be  written  as 

(c«  +  AJB)n«  +  X)  (^fc***  +  »  (2*1) 

h 

where  represents  solvent-solvent,  solvent-surface,  and  all  other  interactions  which  do 
not  e:q>Ucitly  involve  electronic  degrees  of  freedom,  is  the  vacuum  energy  level  of  the 
acoq>tor  km’s  electron  orbital  whidi  is  involved  in  the  charge  transfer  (orbital  **0”),  and  A£7 
is  the  shift  of  the  energy  levri  due  to  the  presence  of  the  fluctuating  solvent.  The  acceptor 
orbital  has  an  occupant  equal  to  either  0  or  1.  The  sum  over  k  in  Eq.  (2.1)  rq>resents  a 
sum  over  the  electronic  energy  states  cfc  in  the  semi-infinite  metal  electrode.  This  summation 
term  clearly  includes  contributions  arising  from  both  the  occupancy  of  the  metal  election 
states  as  well  as  a  piece  which  is  a  ‘Hransfer  term”  giving  rise  to  charge  transfer  between  the 
states  and  the  km  orbital  (c^  and  c  are  creation  and  annihilation  operators  respectively). 

The  Hamiltonian  of  Eq.  (2.1),  with  its  explicit  quantum  operators,  is  not  yet  in  a  useful 
fonn  for  the  purposes  of  computer  simulation.  To  transform  it  to  the  desired  form,  one 
can  follow  the  wink  of  Grimly  [16]  and  Muscat  and  Newns  [17].  First,  the  physically 
reasonable  assumption  is  made  that  the  solvent  adiabaticaUy  follows  the  electronic  degrees 
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of  fireed(»ii  due  to  the  different  timescales.  One  can  then  examine  jtist  the  electronic  part  of 
the  Hamiltonian  for  a  given  solvent  shift  variable  A  i?  treated  adiabatically  : 

Ha  —  (ca  +  Ai?)na  +  5^  (cfcnfc  +  V^c^Cfc  +  VioCfcOa)  •  (2.2) 

h 

By  u  series  of  mathematical  maniptdations  involving  the  resolvent  operator  for  this  Hamil- 
tonian,  it  can  be  shown  that  the  ground  state  electronic  energy  is  given  up  to  a  constant  by 
the  i^ptofximate  form  [16,17] : 

£i(A£)=  ^  +  i(t.  +  Ag  -  v)ton-‘  (V  ~ 

+  ^ln((<,  +  Ag -«/)’  + A”)  ,  (2.3) 

where  e/  is  the  Fermi  level  of  the  metal  and  A  is  a  parameter  which  measures  the  broadening 
of  the  energy  of  the  orbital  a  due  to  the  presence  of  the  metal  riectrons.  The  parameter  A 
is  formally  given  by  the  expression  : 

A(e)  =  ir5:iKa,|»5(€~efc)  .  (2.4) 

k 

Here,  Vak  is  the  matrix  element  coupling  the  orbital  and  a  particular  metal  electron  state  k. 
As  it  is  roughfy  insensitive  to  the  energy  e  for  energies  in  the  range  of  interest,  A  is  typically 
taken  as  a  constant. 

The  net  effect  of  the  preoeeding  manipulations  is  that  classical  dynamics  can  now  be  run 
on  the  Hamiltonian  HmIv  +  £o(  A£7)  which  represents  the  adiabatic  ground  state  behavior  of 
the  Hamiltonian  of  Eq.  (2.1).  This  result  is  obtained  because  one  can  easily  assign  AE  to 
certain  classical  potential  energies  in  a  typical  molecular  dynamics  simulation.  In  partiailar, 
in  the  Hamiltonian  fiw  the  present  simulation  one  can  separate  out  all  the  terms  that  depend 
on  the  Colonmbic  interacticm  between  the  s^dvent  and  the  ion.  One  can  then  write  those 
terms  in  a  form  assuming  a  ferric  ion  (Fe*^)  with  a  positive  charge  Z  (Z  =  3)  which  has  an 
electom  orlntal  of  occupancy  Ua.  ff  the  orbital  is  empty,  then  Ua  =  0,  and  the  Hamiltonian 
effactiviriy  reduces  to  that  describing  the  Fe*^  ion  solvated  near  the  metal  sur&ce.  ff  the 
km  orirital  is  filled,  however,  then  =  1  and  the  Hamiltonian  describes  the  combination 


4 


of  the  Fe^  ion  with  the  negatively  charged  electron,  i.e.,  a  solvated  ferrous  ion  (Fe^'*' )  near 
the  surface.  The  overall  Coloumbic  terms  are  thus  given  by  : 

{Z  -  (2.5) 

Hoe,  Vtjhn  indicates  the  potential  energy  between  the  solvent  and  a  single  positive  charge 
at  the  urn’s  position,  while  indicates  the  potential  energy  between  the  solvent  and  a 

unit  negative  charge  at  the  position  of  the  k>n’s  image.  In  order  to  make  a  correspondence 
between  Eq.  (2.5)  and  the  Anderson-Newns  Hamiltonian  of  Eq.  (2.1)  one  must  therefore 
specify 


(2.6) 

The  Z'dependent  terms  in  Eq.  (2.5)  are  grouped  with  the  rest  of  the  interactions  in  the 
solvent  Hamiltonian  Htohf 

The  end  result  of  the  analysis  is  that  adiabatic  computer  simulations  can  be  performed 
for  an  ion  and  its  image,  with  an  additional  term  in  the  Hamiltonian  rq>re8enting  the 
eJectrooic  degrees  freedom  given  by  Eq.  (2.3).  One  therefore  expects  the  free  energy 
curves  wiU  display  the  approximate  form  [21,19,20] : 

F(Ae)  =  -  AEUm)*  +  Eo(^J5)  (2.7) 

Here,  A  is  the  reorganization  free  energy  obtained  from  a  best  fit  parabola  to  the  diabatic 
Fe^  free  energy  curve,  and  AEn»«n  gives  the  location  of  the  minimum  of  the  diabatic  curve. 

Standard  umbrella  sampling  techniques  (see,  e.g.,  refs.  [5,8]  can  be  applied  to  calculate 
the  adiabatic  free  energy  curve  as  a  function  of  AE  and  find  the  free  energy  barrier  to  the 
charge  transfer  reaction  involving  the  ionic  orbital  a.  Also,  the  average  occupancy  of  the 
accqitor  orbital  a  at  any  point  along  the  free  energy  curve  can  be  shown  to  be  given  by 

[1«»17] 

MAB))  =  1/2  +  (l/»)tan-‘  h  ~  j  (2^) 
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Equation  (2.8)  serves  as  a  useful  summary  of  the  behavior  of  the  Anderson-Newns  Hamil¬ 
tonian.  The  average  occupancy  of  the  orbital  a  on  the  ion,  (na(AE)),  smoothly  goes  from 
0  to  1  as  a  function  of  the  shift  in  its  energy  due  to  the  solvent,  given  by  AE.  If  the  solvent 
shift  takes  on  a  value  such  that  the  energy  level  is  altered  to  coincide  with  the  Fermi  level  of 
the  metal,  the  average  occupancy  is  one-half  (i.e.,  the  electron  transfers).  The  parameter  A 
measures  how  quiddy  the  transition  from  0  to  1  occurs  as  AE  is  varied,  which  is  in  keeping 
with  its  role  as  an  energy  level  broadening. 

B.  Model 

Some  of  the  simulation  details  in  the  present  work  were  presented  in  an  earlier  article 
by  the  present  authors  [11]  but  are  repeated  again  here  in  the  interest  of  completeness. 
The  solvent  consisted  of  671  water  molecules  described  by  a  flexible  SPC  model  [22,23]. 
This  model  was  modified  in  a  fashion  which  is  inconsequential  for  the  classical  case,  but  is 
iiiq>ortant  for  later  generalization  to  the  quantum  case,  and  is  discussed  in  Appendix  A.  A 
tenq>eratiure  of  300  K  was  maintained  during  the  simulation  by  using  one  Nose  thermostat 
[24]  on  the  oxygens,  and  another  on  the  hydrogens. 

The  water  molecules  in  the  simulation  were  situated  at  the  interface  with  a  Pt(lll) 
dectrode  sur&oe.  The  interactions  between  the  water  molecules  and  the  platimun  surface 
were  ^ven  by  the  potential  developed  by  Raghavan  et  al.  [25]  which,  in  tmm,  was  developed 
by  fitting  to  the  water-platinum  atom  potential  given  by  Spohr  and  Heinzinger  [26].  Since 
this  effective  potential  represents  the  Pt(lll)  plane  with  its  fee  lattice  structure,  the  potential 
ediibits  an  hexagonal  symmetry.  Moreover,  this  potential  represents  the  corrugated  Pt 
(111)  sur&ce  without  including  individual  platinum  atom*;  in  the  simulation.  Using  such 
an  effective  potential  is  advantageous  as  one  does  not  have  to  devote  s^nificant  computer 
resources  to  the  integration  of  the  metal  atoms’  trajectories  which  are  unlikely  to  contribute 
significant  effects. 

Hescagonal  periodic  boundary  conditions  were  not  employed  in  the  present  study  as  in 
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the  study  of  Raghavan  et  al.  [25]  Rather,  a  carefully  chosen  rectangular  cell  was  employed 
whose  sides  of  length  3s  and  y/2a  can  inscribe  a  hexagon  of  side  length  a,  shown  in  Fig.  1. 
Periodicity  was  enforced  in  the  x  and  y  directions.  Such  a  simulation  cell  can  accommodate 
the  underlying  hexagonal  symmetry  of  the  (ill)  surface,  as  shown  in  Fig.  2.  For  the 
particular  case  of  the  (111)  face  of  platinum,  the  snoallest  hexagon  that  could  be  formed 
had  a  side  of  length  2.77  A.  A  central  hexagon  with  a  side  five  times  larger,  or  13.86  A,  was 
chosen  which  gave  a  box  length  in  the  x  direction  of  41.6  A  and  in  the  y  direction  of  24.0  A. 
The  water-water  interactions  were  smoothly  cut-ofT  by  half  the  y  box  length,  going  to  zero 
at  12.0  A.  The  smooth  cut-off  b^an  at  a  distance  of  10.3  A. 

While  the  bottom  of  the  periodic  cell  was  bounded  by  the  platinum  surface,  the  top  of 
the  cell  was  bounded  by  a  featureless  slab.  The  parameters  for  the  interaction  of  this  slab 
with  a  water  molecule  were  taken  from  the  ‘‘soft**  water-metal  potential  of  Hautman,  Hallqr, 
and  Rhee  [28],  with  the  image  terms  omitted.  The  hdght  of  this  slab  was  adjusted  to  give 
bulk  water  density  in  the  middle  5  A  of  our  cell,  giving  the  periodic  cell  a  length  in  the  z 
direction  of  22.5  A. 

Raghavan  et  al.  [25]  have  demonstrated  that  their  potential  leads  to  a  rather  interesting 
water  structure  at  the  metal  liquid  interface  for  SPC/E  water  [29].  Since  the  distance 
between  the  platinum  atoms  on  the  (111)  face  2.77  A)  is  conunensurate  with  the  distance 
between  water  molecules  in  ice,  the  waters  apparently  form  an  ice-like  layer  on  top  of  the 
Pt(lll)  face  which  esdiibits  hexagonal  ordering  even  at  room  temperature.  There  is  also  a 
second  layer  of  waters  which  is  hydrogen  bonded  to  the  first  one.  Beyond  the  second  layer, 
however,  the  water  exhibits  essentially  bulk-like  characteristics  (although  there  is  some  hint 
of  a  moderate  nonuniformity  in  density  which  could  be  interpreted  as  corresponding  to  a 
third  partially  ordered  layer).  It  shotdd  be  noted  that  even  though  a  flexible  water  model 
has  been  employed  in  the  present  study  which  is  slightly  different  from  rigid  SPC/E,  it  was 
again  verified  that  the  water  interaction  with  the  platinum  surface  again  gives  rise  to  a 
similaT  water  double  layer  structure. 

The  iron  ion  was  fixed  5.1  A  above  the  center  of  the  platinum  surface  which  places  it  in 
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the  center  of  the  second  layer  of  waters.  This  placement  is  in  accordance  Miith  the  work  of 
Rose  and  Benjamin  [13].  The  interaction  between  the  iron  zind  the  water  was  taken  from 
Kuharsld  et  al.  [8].  It  essentially  consists  of  a  Coloumbic  interaction  between  the  charge  on 
the  ion  with  the  partial  charges  on  the  SPC  water,  along  with  a  repulsive  potenti2d  between 
the  iron  and  the  oxygens.  In  addition  to  this  iron  ion,  an  image  charge  of  the  ion  was 
intoduced  below  the  surface,  as  discussed  in  Appendix  B.  All  iron-water  interactions  were 
cut-off  in  the  same  fashion  as  were  the  water-water  interactions. 


C.  iVee  Energy  Sampling  Method  -  Umbrella  Sampling 

The  free  energy  curve  F(Ae),  is  related  to  the  probability  of  observing  a  certain  value  of 

the  energy  gap  function  AE  =  Ae.  Therefore  one  could  run  simulations  in  equilibrium  and 

via  some  sort  of  binning  arrangement  calculate  the  probability  of  obtaining  various  values 

of  Ae  and  thereby  extract  F{Ae).  However  the  reactive  ET  event  occurs  very  infrequently. 

Consequently  one  would  be  unable  to  effectively  sample  these  important  transitional  regions 

of  the  free  energy  curve  using  a  realistic  amount  of  computer  time.  The  solution  to  this 

problem  is  to  add  a  biasing  potential  to  force  the  system  to  obtain  a  particular  range  of 

values  of  Ae  and  then  remove  the  effect  of  the  potential  when  calculating  F(Ae).  One  can 

see  how  to  do  this  by  simply  rewriting  the  equilibrium  average  of  some  quantity  A  in  the 

system  described  by  a  Hamiltonian  H  (i.e.  {A)h)  In  terms  of  a  different  system  with  the 

added  biasing  potential  H  : 

■pH  A  r  A  f 

(2.9) 


...  /  <Jx  e-^A  s  <ix  /  dx 

\A^g  .  .  ^ww  .  .  ^wr  X 


s  dxe-^«  ~  s  dx  e-^»  '  S  dx 

{A)„  =  (2.10) 

Specializing  to  the  case  A  =  S{Ae  —  AE)  and  H  =  Vf,{AE),  and  using  the  definition  of 
the  free  energy  (up  to  an  overall  constant),  one  obtains 

F(Ae)  =  -keT  ln(5(Ae  -  AE))h  (2.11) 

=  -fcbrin(5(Ac-AE))H+H-Vfc(Ae)-bc  .  (2.12) 
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Hoe  c  is  some  constant  that  depends  on  the  specific  choice  of  the  biasing  potential  By 
using  a  series  of  biasing  potentials  in  overlapping  regions  of  the  variable  Ae,  i.e.  ‘^windows”, 
one  can  match  up  the  results  and  form  a  smooth  curve  of  F(Ae)  for  the  original  system, 
even  in  otherwise  infrequently  sampled  regions  relevant  to  the  ET  process.  [5,8] 

D.  Reactive  Flux  Calculations 

In  the  adiabatic  limit,  the  Marcus  rate  constant  is  essentially  a  classical  transition  state 
theory  expression  obtained  by  using  Ae  as  the  reaction  coordinate.  [1]  In  other  words,  the 
theory  ignores  dynamical  recrossings  so  that  once  the  reaction  coordinate  passes  through  the 
transition  state  towards  the  product  it  is  assumed  that  it  will  never  turn  around  and  come 
back.  A  standard  technique  known  as  the  reactive  flux  correlation  function  [30]  has  been 
developed  to  efficiently  calculate  the  efect  of  these  recrossing  -  an  effect  which  is  expressed 
in  terms  of  a  piefactor  k  to  the  TST  rate  constant.  The  reactive  flux  correlation  function 
expression  for  k  has  its  origin  in  the  classical  expression  for  the  rate  constant  to  go  from  a 
stable  state  A  (reactant)  to  a  stable  state  B  (product)  : 

=  (*x)-‘«(0)«(«(0) -, •)««(,(«))>  .  (2.13) 

Here,  q  is  the  generalized  reaction  coordinate,  q*  denotes  its  transition  state  value,  6b  is  a 
unit  step  function  with  its  origin  at  the  transition  state,  and  denotes  the  molar  fraction 
in  the  reactant  state  (x^  +  =  1).  The  zero  time  limit  of  this  expression  gives  the  TST 

rate  constant  [30]  which  assumes  no  recrossings  and  is  an  upper  boimd  to  ksA  • 

krsT  =  (*x)~^(9(0)5(q(0)  -  g*)0B(q(O)))  .  (2.14) 

Since  the  TST  quantity  is  an  upper  botmd,  one  can  write  the  exact  rate  constant  as 

ksA  =  «  krsT  >  (2.15) 

where  /e  <  1.  The  closer  /c  is  to  unity  the  less  important  are  the  transition  state  recrossings. 
With  (...)«  denoting  here  the  constraint  q  =  q*  for  the  initial  conditions,  one  can  write 
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freedozxL  Tlie  free  energy  as  a  function  of  the  solvent  coordinate  is  defined  up  to  a  constant 
by  setting  the  probability  of  finding  a  certain  value  Ae  of  as  equal  to  exp[— /9F(Ae)]. 
More  explicitly,  one  has  [5,4]: 

F(Ae)  =  -fcBrfn[ydx^(A£?(x)- Ae)e-'5’"(*)j  ,  (3.1) 

where  x  represents  the  mvitidimensional  coordinates  of  the  entire  system. 

Anal<^ous  to  the  discussion  in  our  earlier  paper  [11],  the  reaction  coordinate  can  Ls 
written  as  the  negative  of  the  difference  of  the  solvent  interacting  with  the  3+  ion  and  the 
2+  ion,  since  the  non-ootilombic  interactions  are  identical  in  both  cases.  Using  this  fact  and 
Eq.  (3.1),  it  can  be  shown  [4]  that  the  free  energy  curves  satisfy  : 

Fj+(Ae)  =  F3+(Ae)  +  Ae  (3.2) 

As  a  result  of  the  above  argument,  only  the  2+  diabatic  free  energy  curve  was  calculated 
by  the  techniques  of  umbrella  sampling  [5,8],  and  the  34-  curve  was  generated  via  Eq. 
(3.2).  Parabolic  biasing  potentials  on  the  solvent  coordinate  AF(x)  were  used  to  define  the 
umbrella  sampling  windows,  and  each  window  was  first  equilibrated  for  over  10  picoseconds. 
Then  the  data  were  collected  for  each  window  from  8  ps  trajectories.  A  total  of  62  partially 
overlapping  windows  were  employed  to  generate  the  curver  shown  in  Figure  3. 

Before  computing  the  adiabatic  free  energy  curves,  one  needs  to  use  the  diabatic  curves 
to  find  the  reorganization  energy  A  and  the  location  of  the  minimum  of  the  Fe^  curve, 
Aerntn*  It  has  been  shown  by  Schmickler  [19]  that  if  the  diabatic  free  energy  curves  are 
parabolic  the  adiabatic  free  energy  curves  will  have  the  form  of  a  synunetric  double  well 
provided  the  overpotential  77  is  set  to  zero.  The  latter  quantity  is  defined  by 

»/  =  -  €/  +  Ae^  -  A  .  (3.3) 

The  quantity  Cq  —  e/  is  adjustable  by  simply  varying  the  external  potential  of  the  electrode, 
but  one  must  know  the  other  quantities  in  the  above  expression  in  order  to  affect  the  desired 
adijustment  to  achieve  the  desired  sjmimetric  form.  Assuming  parabolic  behavior,  the  3-f 
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curve  should  s&tisfy  the  equation  y  =  l/4A(Ae  —  Acmm)^  +  and  a  least  squares  fit  to 
the  classical  data  in  Figure  3  gives  A  =  57.3  kcals  and  Ae^m  =  449.7  kcals,  each  correct 
to  about  0.5  kcals.  Thus,  there  will  be  an  error  of  approxixnately  0.7  kcal  in  the  effort  to 
acheive  a  ssero  overpotential. 

B.  Adiabatic  Curves 

FVv  a  given  overpotential  (e.g.  zero  in  the  present  case),  one  can  determine  the  adia¬ 
batic  curves  using  the  Anderson-Newns  Hamiltonian  (i.e.,  with  Eq.  (2.3)  included  in  the 
dynamics)  and  umbrella  sampling.  Following  the  same  techniques  as  described  previously, 
the  adiabatic  curve  was  found  for  A  =  0.05  e.v.  as  shown  in  figure  4.  Twenty  eight  partially 
overlapping  sampling  windows  were  used  to  compute  the  curve.  Once  the  adiabatic  curve 
was  calculated,  the  top  of  the  barrier  was  found  to  be  at  Ae*  =  390.5  kcals  ±  0.5  kcals. 
By  using  this  information,  reactive  flux  techniques  could  then  be  used  to  determine  the 
inqxxrtance  of  dynamical  recrossing  effects  in  the  ET  rate  constant.  [13] 

C.  Becrossiiig  Behavior 

As  discussed  earlier,  the  effects  of  classical  recrossings  of  the  transition  state  are  measured 
a  quantity  k  given  by  Eq.  (2.19).  In  that  equation,  Ae*  indicates  the  transition  state 
value  of  AE  and  the  symbol  (. .  .)^  denotes  an  equilibrium  average  over  all  the  individual 
solvent  momenta.  The  symbol  (. .  denotes  an  equilibrium  average  over  all  the  individual 
solvent  coordinates  with  the  multidimensional  reaction  coordinate  AE  constrained  by  a  £ 
functk>n  to  be  at  Ae*.  In  practice  this  quantity  was  determined  by  running  trajectories  held 
in  the  area  of  the  barrier  top  by  a  constraining  potential  whidi  is  a  function  of  AE.  When 
AE  passed  through  the  barrier  top,  that  solvent  coordinate  configuration  was  saved  and  the 
solvent  momenta  were  reselected  firom  the  appropriate  Boltzmann  distribution.  A  number 
of  such  trajectories  with  reselected  momenta  were  then  run  and  the  nested  averaging  was 
petfeirmed  over  these  trajectories  with  different  momenta  (the  Nose  thermostats  were  turned 


off  for  this  procedure).  More  runs  were  then  performed  for  different  solvent  configurations 
at  the  barrier  top,  so  as  to  arrive  at  a  final  converged  value  of  k. 

Owing  to  the  discrete  timestep  in  the  simulations,  the  values  of  AE  which  were  saved 
were  very  slightly  off  the  exact  chosen  barrier,  whether  to  the  right  or  the  left.  The  net  effect 
of  this  numerical  uncertainty  is  that  a  delta  function  constraint  is  not  predsly  enforced  for 
the  reaction  coordinate  at  the  top  of  the  barrier.  (This  is  allowed  in  the  reactive  flux 
method  provided  the  barrier  is  high  and  the  stable  states  are  well  separated.  [30])  One  can 
then  deduce  foun  Eq.  (2.19)  that  the  zero  time  value  of  k  will  be  »  0,  but  it  rises  quickly 
and  subsequently  foils  off  to  the  final  plateau  value.  The  graph  of  k  versus  time  is  shown  by 
Fig.  5  for  the  present  qrstem.  In  the  calculations,  112  different  starting  configurations  of  the 
solvent  were  used  to  generate  the  graph,  each  one  being  run  with  100  reselected  momenta 
configurations.  The  final  value  for  k  is  0.5769  +/-  0.0004.  This  shows  ^^'at  recrossing  effects 
do  not  substantially  cast  a  shadow  on  the  TST  (Marcus)  estimate  of  the  ET  rate  constant 
whidi  describes  the  rate  over  many  orders  of  magnitude. 

IV.  QUANTUM  EFFECTS 
A.  Path  Integral  Methods 

Since  classical  recrossings  are  found  to  have  little  effect  on  the  overall  rate  constant, 
the  role  of  quantiun  mechanical  solvent  effects  was  next  examined.  An  earlier  study  of  the 
homogeneous  Fe^/Fe^  ET  system  found  tunneling  of  the  water  librational  modes  can  can 
have  a  substantial  impact  on  the  rate  constant,  [9]  although  it  has  been  suggested  [10]  that 
the  latter  effect  may  be  somewhat  enhanced  by  the  particular  model  used  for  water.  Clearly, 
the  issue  of  solvent  quantization  is  an  imp(»tant  one  in  ET  theory,  and  the  present  study 
examines  these  effects  for  the  first  time  for  the  quite  different  case  of  heterogeneous  ET. 

The  quantum  mechanical  simulations  woe  carried  out  using  the  well-established  tech¬ 
nique  of  path  integral  molecular  dynamics  [32,33].  In  order  to  use  such  a  simulation  to 


actually  calculate  an  activated  rate  constant,  path  integral  quanttim  TST  was  \ised.  [18] 
This  theory  provides  a  quantum  medianical  analogue  for  the  classical  TST  activation  free 
energy  by  expressing  the  quantum  rate  constant  k  as 

k  —  K  kQTST 

where  fc  is  a  quantum  correction  factor  of  order  imity  (but  not  necessarily  less  than  unity), 
and  korsT  can  be  written  in  the  approximate  form 

Here,  is  the  difference  in  free  energy  from  the  well  to  the  barrier  of  the  quantum  free 
energy  curve  corresponding  to  the  path  centroid  of  the  reaction  coordinate  AE.  In  the 
present  case,  the  reaction  coordinate  centroid  is  defined  as 

The  quantity  md  in  Eq.  (4.2)  is  the  frequency  of  the  reactant  well  corresponding  to  the 
centrmd  themud  motion  in  that  r^on.  The  transition  state  configuration  in  the  path 
int^ral  QTST  is  given  by  AEd  =  Ae*,  i.e.,  the  maximum  of  the  quantum  free  energy  curve. 
The  adiabatic  free  energy  curves  of  the  centroid  of  AE{r)  will  therefore  be  the  subject  of 
interest  in  the  following  discussion.  The  notation  Ace  will  be  used  to  refer  to  a  particular 
value  of  the  centroid  AEo  when  discussing  the  quantum  simulation;  the  context  of  the  usage 
should  remove  any  potential  confusion  with  the  classical  case. 

In  path  int^jral  sinmlations,  the  isomorphic  quasiparticle  polymer  which  represents  a 
quantum  partide  [33]  is  spread  out  in  spac»,  but  collapses  towards  a  point  partide,  and 
becomes  mote  dassical-like,  in  the  limit  of  high  temperature  or  large  mass.  As  a  consequence, 
only  the  hydrogens  of  the  water  modd  were  discretized;  the  oxygens  are  over  an  order  of 
magnitude  mcne  massive  and  were  treated  classically. 

The  mass  oi  the  quasipartides  in  the  discretized  path  int^ral  is  an  arbitrary  parameter. 
In  practice,  however,  during  molecular  dynamics  simulations  one  would  prefer  larger  vdoc- 
ities  so  that  the  hydrogen  polymers  more  evenly  sample  the  available  phase  space.  This 
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can  be  acoomplished  by  choosing  smaller  quasiparticle  masses,  but  masses  too  small  give 
vdodties  so  large  that  they  pose  difficulties  for  the  numerical  integrator  and  lead  to  non- 
ergodk  saiiq>ling.  A  mass  of  50  a.u.  was  found  to  be  suitable  on  both  fronts.  In  addition, 
one  needs  to  select  a  value  of  the  discretization  parameter  P  (i.e.,  the  number  of  quasipar¬ 
ticles).  [18]  Ideally,  the  value  of  P  selected  would  be  the  smallest  one  which  exhibits  the 
proper  convergence  in  the  interest  of  coiiq>utation  effort.  In  order  to  estimate  the  value 
of  P  which  gives  such  convergence,  a  Monte  Carlo  study  was  performed  by  Lobaugh  [34]. 
Hiis  study  was  along  the  lines  of  one  suggested  by  Kuharsld  and  Rossky  [35].  It  consisted 
of  two  solvent  water  molecules  which  interacted  using  the  flexible  SPC  potential  and  an 
additional  quadratic  constraining  potential  between  their  centers  of  mass.  This  potential 
had  a  minimum  at  2.85  Aand  a  frequency  u>  =  26  ps~^,  creating  an  environment  similar  to 
the  bulk  solvent.  A  total  of  10^  trial  Monte  Carlo  moves  were  carried  out  and  a  number  of 
different  properties  were  examined  as  a  function  of  P,  including  radial  distribution  functions 
and  the  bond  length  average  (vob).  The  value  P  was  selected  as  exhibiting  sufficient 
convergence  in  these  properties;  for  example,  a  graph  of  (^ob)  vs.  P  is  shown  in  Fig.  6. 

The  application  of  the  path  integral  discretization  scheme  to  a  classical  model  for  water 
involves  the  assumption  that  the  classically  optimized  water  model  will  still  provide  a  good 
representation  for  actual  water  in  the  quantum  limit.  One  drawback  of  this  approach  is 
that  the  parameters  of  the  water  model  urare  developed  classically  to  model  actual  water. 
Fortunately,  at  room  temperature,  the  classically  accessible  r^ions  provide  the  dominant 
contribution  to  the  physical  properties.  Tedmically,  however,  the  model  should  be  reparam- 
eterized  in  the  quantum  limit.  Antother  drawback,  which  has  also  not  yet  been  fully  investi¬ 
gated  due  to  the  computational  cost,  is  that  the  quantum  version  of  the  water  may  explore 
n^ions  of  the  potential  energy  surface  that  are  clasucally  inaccessible.  Thus,  the  functional 
form  of  the  potential  in  these  reports  may  not  be  correctly  parameterized  to  reproduce  the 
physical  properties  of  water.  For  ocample,  it  was  found  that  the  discretized  version  of  the 
rbisical  model  dissociates  at  infrequent  intervals  as  the  potential  is  not  bounded  along  the 
0-H  bonds.  This  resulted  from  occasional  unphysical  tunneling  of  an  hydrogen  to  another 
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nearby  water  molecule  (roughly  1  dissociation  occurred  for  every  500  water  molecules  in  a 
picoieoond  of  simulation).  To  prevent  this  behavior  from  occurring,  the  potential  was  modi¬ 
fied  by  the  addition  of  a  hard  repulsive  wall  along  the  0-H  bond.  As  discussed  in  Appendix 
A,  this  wall  was  adjusted  to  prevent  the  dissociations  without  affecting  the  thermodynamic 
properties  of  either  the  quantum  or  classical  versions  of  the  water. 

The  numerical  techniques  used  for  the  quantum  case  were  essentially  identical  to  the 
rlassiral  case.  Again  two  Nos^  thermostats  were  used;  one  on  the  oxygens  and  one  on 
all  of  the  hydrogen  quasipartides.  The  sol^^t  coordinate  consists  of  the  same  Coulombic 
interactions,  with  the  exception  that  eadi  hydrogen  quasiparticle  acts  to  create  a  potential 
at  the  ion  site  as  would  a  charge  of  magnitude  -f-1/.P.  The  same  value  of  —  e/  was  used 
in  order  to  &irly  compare  the  quantum  adiabatic  curves  with  the  classical  ones  (i.e.  to  be 
aUe  to  clearly  discern  quantum  solvent  effects  for  a  single  underlying  Hamiltonian). 

B.  Quantum  Results 

The  results  frmn  the  quantum  umbrella  sanq>ling,  which  contains  contributions  from  9 
sanqding  windows,  are  shown  in  Fig.  7.  It  should  be  noted  that  with  P  =  25  the  simulation 
requires  roughly  11  times  more  computer  time  than  the  classical  case  to  conq>lete  a  given  run. 
In  fiict,  the  present  simulation  seems  to  be  the  most  ambitious  simulation  of  quantised  water 
attenq)ted  to  date.  All  modes  the  fiesdble  model  have  been  quantised  and  a  relativdy 
large  number  of  water  molecules  were  used  (671).  The  simulation  required  the  dedicated  use 
of  five  high  end  computo  workstations  (e.g.  HP/ Apollo  735)  for  4  months  along  adth  200 
C9PU  hours  on  a  Cray  C90.  Of  course,  in  the  future  such  computations  will  become  mote 
ooauDonplaoe  as  comqmters  increase  in  speed. 

As  can  be  seen  in  Fig.  7,  the  quantisation  of  the  solvent  model  has  a  significant  effect  on 
the  ficee  energy  curves.  The  same  value  for  —  Ca  was  used  as  in  the  classical  case,  but  now 
the  carves  are  no  longer  symmetric.  The  overpotential  in  the  quantum  case  is  10.3  kcals, 
creating  free  energy  curves  which  are  in  marked  contrast  to  the  symmetric  classical  case.  This 
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bduKvior  xeiultc  becftiue  the  minimum  reaction  coordinate  value  is  shifted  from  449.7  kcals 
to  464.3  heals  in  the  quantum  case,  while  A  has  shifted  from  57.3  kcals  to  61.6  kcab.  The 
quantum  barrier  hei|;ht  from  the  2+  well  is  8.9  kcals,  as  opposed  to  11.2  kcals  in  the 
daancal  case.  For  the  ~  3+  wdl,  the  quantum  barrier  is  19.1  kcals  as  opposed  to  the  classical 
12.7  kcals.  The  results  illustrate  that  the  same  electrode  overpotential  will  not  necessarily 
result  in  wmilay  adiabatic  free  energy  curves  in  the  classical  and  quantum  mechanical  models. 
In  fact,  the  thermodynamic  driving  force  is  substantially  different  for  the  two  cases,  reflecting 
the  differences  in  sedvation  of  the  ion  by  the  classical  and  quantum  water  modds.  These 
differences  in  free  energies  would  have  quite  large  effects  on  the  rates.  This  study  therefore 
strongly  suggests  the  inq>ortanoe  of  treating  quantum  solvent  effects  within  the  quantum 
TST  eaqwession  (i.e.  a  factor  of  60  for  the  oxidation  of  Fe^  and  of  4  x  10^  for  the  reduction 
of  in  studies  of  heterogeneous  electron  transfer  rates  for  aqueous  solutions. 

V.  CONCLUDING  REMARKS 

The  present  studies  have  examined  the  role  of  both  a  classical  and  quantum  mediani- 
cal  water  solvent  modd  with  regards  to  electron  transfer  across  the  water-metal  inter&oe. 
Vl^thin  a  Marcus  theory  framework,  solvent  dynamical  effects  have  been  shown  to  be  of 
rdativdy  minor  importance.  A  comparison  of  the  quantum  adiabatic  free  energy  curves 
with  the  corresponding  classical  curves  has  shown,  oa  the  other  hand,  a  considerable  impact 
from  quantum  effects.  Since  a  quantum  solvent  modd  more  accuratdy  describes  the  actual 
idiysical  qrstem,  thw  result  suggests  that  a  classical  simulation  may  not  be  adequate  for  the 
stwly  of  heterogeneous  ET. 

A  comparison  of  the  present  results  with  those  of  Bader  et  al.  [9]  on  the  quantum  solvent 
effects  in  the  homogeneous  Fe^/Fe*^  ET  system  is  in  order.  In  the  latter  study,  the  path 
integral  method  was  used  to  quantize  the  rigid  SPC  modd  for  water.  In  that  case,  the 
important  Hbrational  motions  of  the  water  dipoles  were  quantized  and  therefore  allowed  to 
tumid  in  the  symmetric  Fe^/Fe^  ET  process.  Indeed,  these  librational  tunnding  effects 
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wm  estimated  to  be  significant  (a  &ctor  of  60),  although  Song  and  Marcus  [10]  suggest  this 
&ctor  ivould  be  revised  downward  by  considering  a  more  realistic  solvent  spectral  density. 
In  the  present  simnlations  of  the  single  ion,  heterogeneous  ET  process,  all  of  the  modes  of 
the  flexible  SPC  model  have  been  quantized.  Themfote  it  is  more  difficult  to  specifically 
attribute  the  quantum  effects  to  librational  tunneling  (although  they  are  surdy  present), 
hfoeeover,  the  heterogeneous  Fc^'^/Fe^  EST  system  is  quite  different  from  the  homogeneous 
case,  since  the  symmetry  of  the  homogeneous  system  is  broken  by  replacing  one  of  the  iron 
ions  with  a  metal  electrode.  In  the  heterc^eneous  case,  the  differences  in  the  classical  versus 
quantum  km  solvation  may  be  the  dominant  effect  whidh,  in  fact,  leads  to  changes  in  the 
thermodynamic  driving  force  for  the  reaction  (c.£.  Fig.  7).  Sudi  effects  are  not  possible  in 
the  symmetric  homogeneous  Fe^^/Fet^  ET  qrstem. 

Fbiture  researdi  must  include  a  thorough  and  costly  examination  of  the  quantum  water 
modd  to  insure  its  accuracy.  It  is  unclear  whether  the  relativdy  straightforward  exten- 
fton  of  a  classical  model  to  the  quantum  r<^;ime  via  a  path  integral  discretization  scheme 
provides  a  suffidently  accurate  representation  o£  physical  water  without  reparameterization 
of  the  modeL  Of  partiduar  concern  are  the  r^ons  of  the  potential  whidi  are  classically 
inaccessible.  Sudi  r^pons  will  undou'jtedly  need  at  least  minor  modification.  Additional 
researdt  must  abo  indude  an  examination  of  the  role  of  solvent  dectronic  polarizability,  as 
wdfi  as  a  serious  attempt  to  develop  an  effective  classical  water  modd  which  can  capture 
the  quantum  solvation  effects.  [36]  The  inqiortant  issues  of  the  water/surface  interaction 
and  the  electronic  melange  term  A  (cil  Eq.  2.4)  clearfy^  require  some  serious  theoretical 
attention,  as  does  the  problem  of  finite  icm  concentrations.  These  and  other  issues  will  surdy 
make  the  thecnetkal  study  and  simulation  of  heterogeneous  ET  processes  a  challei^ing  and 
stimulating  endeavor. 
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APPENDIX  At  THE  HARD  WALL  MODIFICATION  TO  THE  FLEXIBLE  SPC 
WATER  MODEL  IN  THE  QUANTUM  LIMIT 

In  this  Appendix,  a  hard  wall  modification  to  the  flexible  SPC  water  model  is  presented. 
As  discussed  earlier,  it  was  found  that  the  discretized  path  integral  version  of  the  flexible 
SPC  water  model  dissociates  at  infirequent  intervals  since  the  potential  is  not  bounded  along 
the  0-H  bonds.  This  behavior  results  from  occasional  unphysical  tunneling  of  an  hydrogen 
from  a  water  molecule  to  an  oxygen  site  on  a  nearby  water  molecule.  To  prevent  this 
situation  from  occurring,  the  potential  was  modified  by  the  addition  of  a  hard  repulsive  wall 
along  the  O-H  bond,  given  by  the  fundonal  form 

“  (rojr.  -  rs)»  +  (roft  -  r*)“  ’ 

vdieKe  oa  ==  15.13  (kcai/mole)(A)*’  and  ra  =  2.6  where  roHm  is  the  equilibrium  OH 

bond  length  distance.  These  parameters  were  adjusted  to  prevent  the  dissociations  without 
affecting  the  thermodynamic  properties  of  either  the  quantum  or  classical  versions  of  the 
water. 

Several  test  rimnlations  were  carried  out  with  P  =  5  in  a  three  dimensional  periodic 
c|ibe  with  512  water  molecules.  It  was  found  that  for  ra  =  2.8  tob^^  or  larger,  the  waters 
sHB  dissociated.  As  tunneling  effects  will  surely  increase  with  larger  P,  this  behavior  wiU 


19 


ckariy  be  present  fat  P  =  25.  For  a  value  of  ri,  =  2.6  roH.c«  tbe  dissociations  ceased.  It  was 
also  e3q;>]icitly  verified  that  this  value  of  prevents  the  dissocations  with  P  =  25  and  does 
not  a&ct  the  known  properties  of  SPC  water.  In  addition,  the  average  total  intramolecular 
potential  for  the  P  =  5  case  was  examined  as  a  function  of  and  was  found  to  have 
converged  by  ri,  =  2.6  This  behavior  is  shown  in  Fig.  8. 


APPENDIX  Bt  IMAGES  AT  THE  DIELECTRIC-CONDUCTOR  INTERFACE  IN 

COMPUTER  SIMULATIONS 


In  this  Appendix,  the  incorporation  of  image  charges  into  molecular  dynamics  simu¬ 
lations  at  a  dielectric-conductor  inter&oe  is  examined.  This  analysis  is  performed  in  such 
a  way  as  to  make  the  appropriate  connections  with  a  didectric  continuum  model  and  a 
perfectly  flat  conducting  surface. 

ff  an  ion  is  in  solution  near  a  conducting  surface,  one  would  expect  an  image  term  to  be 
present  in  the  Hamiltonian  for  the  system.  In  feet,  it  b  a  classic  problem  to  have  a  charge 
q  imbedded  in  a  dielectric  medium  of  dielectric  constant  Ci,  nc  ur  the  interface  with  another 
medium  ci  dielectric  constant  £2  (the  conducting  case  corresponds  to  the  limit  of  £2  —*  00). 
If  the  half  space  z  >  0  is  filled  with  dielectric  £1,  while  z  <  0  is  filled  with  dielectric  £2, 
Jadcaon  [37]  solves  for  the  potential  4  when  a  charge  q  is  imbedded  in  medium  1  a  distance 
d  feom  the  interface.  Fbr  the  region  z  >  0  the  solution  is  : 


$  = 


-2.  +  _2L 


(Bl) 


where  R  measures  the  distance  from  an  observer  to  the  charge  9,  and  R!  measures  the 
distance  from  an  observer  to  an  image  charge  whidi  is  located  a  distance  d  behind  the 
interface  (and  thus  2d  from  the  actual  charge  g).  The  image  charge  has  a  magnitude  given 
by : 


9^  = 


£2  —  £1 

- g  . 

€a  +  Cl 


(B2) 
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One  poaaUe  interpretation  for  the  potential  given  in  Eq.  (Bl)  is  to  regard  the  expression 
as  representing  two  actual  diarges  q  and  qf  imbedded  in  a  mediiun  of  dielectric  constant  Ci 
whidi  fills  all  of  space,  keeping  in  mind  that  one  can  only  use  the  equation  for  the  z  >  0 
region*  In  the  z  <  0  region  the  potential  is  gi^^  by 

*  =  .  (B3) 

where  R  measures  the  distance  from  the  actual  charge  g,  and  q"  is  given  by 


9"  = 


2e3 

- g 

+  €i 


(B4) 


Fbr  the  case  under  study  in  the  present  work,  medium  2  is  a  conductor,  so  one  is  interested 
in  the  field  in  the  z  >  0  r^on  when  £3  -*  00.  Equation  (B2)  clearly  shows  that  q'  =  — q, 
and  the  potential  from  Eq.  (Bl)  gives  : 


While  Eq.  (B5)  is  the  correct  expression  for  the  field  in  dielectric  1,  it  is  nevertheless 
not  yet  in  a  form  useful  for  direct  incorporation  into  a  molecular  dynamics  simulation.  The 
reason  for  this  lies  in  the  maumer  in  which  such  a  simulation  takes  into  account  the  dielectric 
cmistant.  No  explicit  dielectric  constant  is  ever  introduced  into  the  intermolecular  forces. 
Rather  the  intermolecular  forces  allow  for  an  orientational  polarizability  (or  other  forms  of 
pdaxfoability  as  well,  depending  on  the  model)  which  responds  to  an  applied  electric  field 
in  such  a  manner  as  to  behave  as  a  dielectric  of  constant  For  SPC  water  at  standard 
tenyerature  and  pressure  Ci  has  been  computed  to  be  65  +/-  9  [8]  while  the  experimental 
water  value  is  80.  It  is  seen  that  even  though  the  expression  for  the  potential  from  a  cha^  in 
bulk  dielectric  may  be  q/eiR,  we  would  carry  out  this  sinnilation  on  a  computer  by  sinq>ly 
aUofwing  the  solvent  to  rdax  around  the  charge  without  any  direct  introduction  of  a 
proper^  omstructed  solvent  model  will  already  take  into  account  its  dielectric  nature.  It 
wfil  arise  of  its  own  accord.  Correspondingly,  if  one  wished  to  carry  out  a  simulation  which 
wodd  give  the  potential  shown  in  Eq.  (B5),  one  cannot  simply  introduce  an  image  charge 
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and  an  explicit  didectric  constant.  Rather,  we  wotild  have  to  do  a  bulk  simulation  where 
the  solvent  is  extended  to  the  2  <  0  region  and  surrounds  a  charge  9'.  This  is  undesirable 
for  several  reasons.  First  of  all,  it  eliininates  not  only  the  water  double  layer  at  the  platinum 
sur&ce,  but  also  the  platinum-water  interaction  entirely.  Secondly,  it  entails  a  much  larger 
simulation  which  is  therefore  substantially  slower  (computer  time  to  perform  the  simulation 
is  roughly  proportional  to  the  square  of  the  number  of  atoms  and  hence  grows  quickly).  And 
finally,  it  wouM  obscure  any  effects  in  the  solvent  polarization  fluctuations  that  might  arise 
due  to  the  solvent’s  restiction  to  the  half-space  z  >  0. 

There  is  a  way  around  this  problem,  however.  One  can  construct  a  simulation  other 
than  the  bulk  one  just  described  which  preserves  both  the  interfadal  water  interaction,  as 
well  as  the  potential  given  in  Eq.  (B5).  This  can  be  done  by  retaining  the  actual  charge  q 
which  is  fixed  a  distance  d  firom  the  interface,  and  introducing  an  additional  image  charge 
g'"  a  distance  d  behind  the  interface  (i.e.  at  the  same  position  as  the  earlier  —  ~q  charge) 
whidi  is  considered  as  being  in  vacuum.  In  other  words,  one  has  a  solvent  with  a  charge  q 
inside  of  it,  with  which  it  interacts  Coloumbically  (and  perhaps  with  other  forces  as  well, 
such  as  Lennard-Jones,  depending  on  the  modd).  In  addition,  the  solvent  interacts  with 
some  effective  potential  that  represents  a  surface  at  z  =  0.  But  there  is  also  an  isolated 
charge  g"'  a  distance  d  bdow  the  z  =  0  plane,  which  does  not  interact  with  the  surface, 
but  does  Coloumbically  (and  only  Coloumbically)  interact  with  the  solvent.  This  charge  g"' 
can  then  be  adjusted  so  that  in  a  continuum  model  it  will  give  the  proper  potential  which 
matdies  with  Eq.  (B5).  Such  a  simulation,  with  an  isolated  charge  behind  the  interface 
interacting  Coloumbically  with  the  solvent,  poses  no  computational  difficulties. 

If  one  rdks  on  the  superposition  of  the  image  and  actual  charges,  the  first  thing  needed  is 
to  write  down  the  potential  in  a  continuum  naodel  given  by  the  actual  charge  in  the  dielectric 
near  another  didectric  of  constant  £3  =  1  (the  vacuum).  The  earlier  Eqs.  (Bl)  and  (B2) 
clearly  show  that  the  potential  in  this  case  is  given  by 

g. _ g_  ,  1  Cl  -1 
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Now  one  can  consider  the  separate  and  distinct  case  of  the  image  charge  cl"  which  is 
located  behind  the  interface  in  the  z  <  0  region.  One  wants  to  find  the  potential  in  a 
continuum  model  that  this  charge  in  vacuum  induces  in  the  z  >  0  region.  To  do  so,  one 
needs  to  adapt  Eqs.  (B3)  and  (B4),  originally  developed  for  z  <  0.  The  rezison  for  this  is 
that  what  matters  is  not  really  whether  z  is  greater  or  less  than  0,  but  rather  whether  one 
is  looking  for  the  potential  in  the  region  that  contains,  or  does  not  contain,  the  charge  of 
interest.  Therefore,  in  a  continuum  model,  the  potential  in  dielectric  1  due  to  the  charge  q"* 
in  vacuum  behind  the  interface  is  given  by 


2c, 


•r . 


etH'  1  +  e, 

Summing  £)qs.  (B6)  and  (B7)  and  setting  the  result  equal  to  Eq.  (BS),  one  finds 


(B7) 


• 


(B8) 


Surptisingiy  the  needed  charge  is  simply  the  negative  of  the  actual  charge,  independent  of 
the  dielectric  constant  c,! 

To  summarize  what  has  been  done :  if  one  is  performing  a  molecular  dynamics  simulation 
with  a  charge  q  in  a  solvent  a  distance  d  above  a  conductor,  one  should  add  an  image  charge 
— q  a  distance  d  below  the  conductor  which  interacts  with  the  solvent.  The  purpose  of  this 
IHOoeduze  is  to  make  rigorous  contact  with  the  expected  results  of  a  dielectric  continuum 
theory.  In  practice,  however,  for  the  simulations  detailed  in  the  present  work  the  ion  is  5.1 
A  from  the  sur&oe.  At  such  a  distance  the  inclusion  of  the  image  charge  with  the  water 
solvent  (e  »  80)  has  only  a  very  minor  effect. 
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FIGURES 


FIG.  1.  Rectangular  periodic  cell  whidi  permits  hexagonal  ordering 

FIG.  2.  Saiiq)le  two  dimensional  tiling  <£  space  with  rectangular  periodic  cell 

FIG.  3.  Solid  lines  fpves  the  diabatic  free  energy  curves  for  the  iron  ion  5.1  A  from  the  Pt 
tuiftoe.  The  upper  solid  line  is  for  the  Fe^'^  ion^  and  the  lower  solid  line  is  for  the  Fe^  ion. 

FIG.  4.  Classical  Adiabatic  free  energy  curves  for  ferrous-ferric  charge  transfer  at  the  wa- 
ter-Pt(lll)  interface.  Here  A  =  0.05  e.v.  Eq.  (2.3)  was  used  for  the  electronic  energy  contribution, 
an  sgiproadznate  expression  to  the  electronic  portion  of  the  Anderson-Newns  Hamiltonian. 

FIG.  5.  K  vs.  time  (eq.  (2.19))  for  ferrous-fmic  charge  transfer  at  the  water-Pt(lll)  interface. 
The  final  value  for  k  is  0.5769  +/•  0.0004. 

FIG.  6.  {tob)  vb.  P  in  the  monte  carlo  path  integral  water  simulation  of  Lobaugh.  P  =  25 
was  sdected  as  exhiUting  converged  beharior. 

FIG.  7.  Adiabatic  firee  energy  curve  for  ferrous-ferric  charge  trcmsfer  at  the  water^Pt(lll) 
intetfime  using  a  quantum  solvent  is  shown  by  the  s<did  line.  Eq.  (2.3)  was  tised  for  the  dectronic 
energy  contribution,  an  approximate  egression  to  a  portion  of  the  Anderson-Newns  Hamiltonian. 
The  datiied  line  shows  a  best  fit  to  eq.  (2.7).  The  dot-datii  line  shows  the  previous  clasacal  case. 

FIG.  8.  <  V^ttra  >  vs.  rh  in  nudecular  dynamics  study  for  P  =  5. 
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